Density oscillations in trapped dipolar condensates 
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We investigated the ground state wave function and free expansion of a trapped dipolar conden- 
sate. We find that dipolar interaction may induce both biconcave and dumbbell density profiles 
in, respectively, the pancake- and cigar-shaped traps. On the parameter plane of the interaction 
strengths, the density oscillation occurs only when the interaction parameters fall into certain iso- 
lated areas. The relation between the positions of these areas and the trap geometry is explored. 
By studying the free expansion of the condensate with density oscillation, we show that the density 
oscillation is detectable from the time-of- flight image. 

PACS numbers: 03.75.Hh, 03.75.Kk, 67.85. Jk 



I. INTRODUCTION 

Bose-Einstein condensates with dipole-dipole interac- 
tion have drawn significant interests over the past few 
years [l|-l3|. In particular, after the experimental obser- 
vations of dipolar effects in Cr condensates [4] and in 
spin-1 Rb [5] condensates, the study of dipolar quantum 
gases is at the forefront of both theoretical and experi- 
mental research. Armed with the technique of Feshbach 
resonance which tunes the short-range 5- wave scattering 
length, experimentalists have investigated the stability 
and the anisotropic collapse of the Cr condensate [6|, Q • 
In addition, dipolar effects have also been observed in 
condensates of K [8] and Li [9] atoms which only pos- 
sess small magnetic dipole moments. Owing to the large 
electric dipole moment of the polar molecules, recent 
success in making high phase-space-density gas of po- 
lar molecules [10, 11] provides an even more promising 
platform for the study the dipolar quantum many-body 
system. 

The long-range and anisotropic characters of the dipo- 
lar interaction have profound impact on the properties of 
quantum gases. For example, when loaded into a optical 
lattice, the long-range repulsive dipole-dipole interaction 
breaks the translational symmetry of the lattice such that 
the density wave phases can be formed [121, Il3| . In the 
presence of vortice s ^141 or a small perturbation in the 
trapping potential [15|, the radial density profile of the 
condensate was found to develop oscillations. Of partic- 
ular interest is that, even for a singly trapped dipolar 
condensate, a biconcave density profile can be found in a 
pancake-shape trap [16]. In this case, the maximum den- 
sity occurs at the periphery of the condensate, resembling 
the shape of a red blood cell. Wilson et al have proposed 
to detect the biconcave density profile by tuning the con- 
tact interaction to induce local collapse |17| . 

In the present paper, we systematically investigate the 
ground state structure and the free expansion dynam- 
ics of a harmonically trapped dipolar condensate. For a 
given trap geometry, we map out the stability diagrams 
on the parameter plane of the interaction strengths. In 



addition to the biconcave dipolar condensates in pancake- 
shaped traps, a dumbbell density profile can be induced 
in cigar-shaped trap when the dipolar interaction is prop- 
erly tuned. For simplicity, both biconcave and dumbbell 
density profiles are referred to as density oscillations in 
the condensates. On stability diagram, the density oscil- 
lation occurs only when the interaction strengths fall into 
certain isolated areas, which is called as density oscilla- 
tion island (DOI). We also find the relation between the 
position of DOI and the trap geometry by fitting the nu- 
merical results. Finally, by studying the free expansion of 
the condensate with density oscillation, we show that the 
density oscillation can be identified in the time-of-flight 
image. 

We will also consider dipolar condensate trapped in 
box potentials. We show that the rigid walls of box po- 
tentials in general induce more pronounced density oscil- 
lations over a much larger region of the parameter space 
than harmonic potentials. 

The remainder of this paper is organized as follows. In 
Secini we briefly outline our model. In Sec. lIIIi we intro- 
duce the numerical method used to compute the dipolar 
mean-field. The numerical results on both the ground 
state wave function of harmonically trapped dipolar con- 
densates and the corresponding free expansion dynamics 
are presented in Sec. IIVI In Sec. |Vl we discuss the en- 
hanced density oscillation in a box potential. Finally, we 
conclude in Sec. IVll 



II. FORMULATION 

We consider a trapped condensate of N bosons with 
dipole moments dz. Particles interact with each other 
via contact interaction Vq and dipole-dipole interaction 
Vd. The total interaction potential takes the form 



V{r) = coS{r) + fd' 



1 - 3 cos2 



(1) 



where cq = A7:h?as/m with m being the mass of the par- 
ticle and as being the 5-wave scattering length tunable 
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FIG. 1: The stability diagram for A = 1. 



through Feshbach resonance, is the polar angle of r, and 
/ G [— ^,1] is a parameter continuously tunable via an 
orienting field which rotates rapidly around z-axis |19j . 
The inclusion of the factor / makes the dipolar interac- 
tion between polarized dipoles tunable. For positive /, 
the dipolar force is repulsive (attractive) if two dipole 
moments aligned side by side (head to tail). However, 
the opposite becomes true for negative /. Moreover, the 
trapping potential is assumed to be axially symmetric 



f/ho(r) 



-muj 



^x'+y'^ + X^z^ 



where uo^ is the radial frequency and A is trap aspect ratio 
such that the axially frequency of the trap is uOz = AcJx • 
The condensate wave function satisfies the non-local 
Gross-Pitaevskii equation, which, in dimensionless form, 
becomes 
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where g = 4:7TNa/a± and D = f d? / {huj ^a\) are the 
dimensionless parameters, characterizing the strengths of 
the contact and dipolar interactions, respecti vely. The 
dimensionless units adopted here are a^ = \/h/{muj±) 
for length, huj± for energy, ujj_ for time, and y^W/a^ for 
wave function. In Eq. (|2]), the control parameters of the 
system are reduced to A, g, and D. We shall show how 
the condensate stability and the structure of the wave 
function depend on these parameters. 



III. NUMERICAL METHOD 

Following the standard procedure, the numerical solu- 
tion of the ground state wave function can be obtained 
by evolving Eq. ([2]) in imaginary time. Compared to the 
Gross-Pitaevskii equation with only contact interaction, 
the difficulty of solving Eq. (|2]) lies at the calculating the 



mean-field potential of the dipole-dipole interaction 

^d(r) = j dv'Vd(v-v')n{v') 



T-^ 



^d(k)n(k) 



(3) 



where T ^ denotes the inverse Fourier transform, n(r) = 
|'^(r)p is the density of the condensate with n(k) being 
its Fourier transform, and Vd(k) = ^ fcP {3 cos^ 6^ — 1) 
is the Fourier transform of the dipolar interaction poten- 
tial with Ok being the polar angle of k. From the second 
line of Eq. (J3j), it becomes apparent that ^d can be effi- 
ciently calculated in momentum space using fast Fourier 
transform (EFT) [2]. Moreover, taking advantage of the 
cylindric symmetry of the system, Ronen et al. [18] fur- 
ther reduce the 2D Fourier transform in (x, y) plane into 
an ID Hankel transform of order for the ground state 
wave function by integrating over the azimuthal variable. 
This procedure greatly simplifies the 3D calculation to a 
2D one in the p-z plane. 

The fact that FFT implicitly treats the system as a 
periodic one makes our system a 3D periodic lattice of 
condensates. Therefore, the dipolar interactions between 
the artificial periodic copies of condensates introduce er- 
ror to the mean- field potential ^di'^)- To solve this prob- 
lem, Ronen et al. [18] suggest to truncate the dipolar 
interaction such that Vd is nonzero only within a sphere 
of radius i?, i.e.. 
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r<R, 

otherwise. 



(4) 



This truncation of the dipolar interaction has no phys- 
ical consequence as long as R is larger than the extend of 
the actual condensate, but it eliminates the interaction 
between artificial copies of condensates. This procedure 
works well in nearly spherical traps. However, for cigar- 
or pancake-shaped trapping potential, one cannot find a 
sphere which only contains a single condensate. Here, a 
natural generalization is to introduce cylindrical trunca- 
tion on the dipolar interaction potential 
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The Fourier transform of V^ ^ (r) takes the form 
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where kn 



the zeroth order Bessel function. The first two lines at 
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FIG. 2: The contact interaction strength dependence of the upper critical dipolar interaction D\ for a condensate in pancake- 
shaped traps. The shaded areas mark the regions of parameter space where the biconcave condensates are found. The broken 
arrows denote the direction to which the DOIs move as A is increased. The inset of (b) shows the typical isodensity surface of 
the condensate in DOT. The inset of (c) reveal the detailed structure of the first two DOIs for A = 13. 



the r.h.s. of Eq. (J6]) are previously obtained in Ref. jlSj , 
and the remaining integral term can be easily worked out 
using numerical integration. With a Gaussian wave func- 
tion, we have confirmed that the precision of the dipolar 
interaction energy is improved significantly by introduc- 
ing cylindrical truncation on Vd- 

Before we conclude this section, let us briefly summa- 
rize the numerical procedure employed to find the ground 
state wave function of a dipolar condensate. For a given 
set of control parameters (A, ^, I)), we first determine the 
proper values of R and Z. Then, using Eq. (J6]), we nu- 
merically compute the Fourier transform of the modified 
dipolar interaction potential Eq. (J6j). Finally, following 
the standard procedure for imaginary time evolution, we 
obtain the ground state wave function. 



IV. DIPOLAR CONDENSATES IN A 
HARMONIC TRAP 



In this section, we shall first map out the stability di- 
agram of a trapped dipolar condensate in the parame- 
ter space {X^g^D). The value of the trap aspect ratio 
A ranges from 0.1 to 20, which covers both highly pro- 
late and oblate trap geometries. For given parameters A 
and D, increasing the repulsive contact interaction will 
stabilize the condensate, we therefore manually set an 
upper limit on g such that g < 600. The lower limit 
of ^, on the other hand, is introduced naturally as the 
system always becomes unstable by tuning the contact 
interaction to attractive with sufficiently large strength. 
Similarly, the lower and upper limits of the dipolar in- 
teraction strength are also determined by the stability of 
the condensate. In the second part of this section, we 
turn to study the free expansion dynamics of an initially 
trapped dipolar condensate, aiming at finding signatures 
of the density oscillation in time-of- flight image. 




FIG. 3: The effective potential t/eff (p, z) of the condensate in 
the presence of density oscillation for (A,^, D) = (8, —10, 20) 
(a) and (0.3,200,-49) (b). 



A. Stability diagram 

In previous numerical studies [sl, [IgI, [22| , the stability 
diagram is plotted in X-D parameter plane for a given 
contact interaction strength g (which is usually taken to 
be zero). Here, for a fixed trap aspect ratio A, we map out 
the stability diagram in g-D parameter plane. As we shall 
show below, this scheme makes it more convenient for us 
to understand the behavior of the density oscillation in 
a trapped dipolar condensate. 

Figure [1] plots the stability diagram in g-D parameter 
plane of a dipolar condensate in a spherical trap, which 
also represents the basic structure of the stability in other 
trap geometries. Since the dipolar force is partially at- 
tractive, a condensate always becomes unstable as one 
increases the absolute value of D. In fact, for given g and 
A, there exist two critical dipolar interaction strengths. 



^(+)/ 



i(-) 
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is stable only if D^^ ^ 



(g) < D < D^^\g). Alternatively, 



if D is fixed, the instability can also be induced by de- 
creasing the scattering length (even to negative), similar 
to what have been done experimentally [3, [21| . 

Moreover, when g (> 0) is so large that the system 
falls into the Thomas- Fermi regime, the variational cal- 
culation shows that the stability boundary for given A 




FIG. 4: The trap aspect ratio dependence of the position 
and the width (represented by the length of the error bar) 
of the DOIs. The hnes are plotted using Eq. (|7| with fitting 
parameters listed in Tab. [H 
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0.2312 


125.6 



TABLE I: The values of fitting parameters for the positions 
of DOIs. 



only depends on the ratio of D and g [22(. This prop- 
erties is clearly reflected in Fig. [TJ as for g > 100, the 
critical D values linearly depend on g with slopes 0.255 
and —0.151 for D^^ and D^^ , respectively. The reason 
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that the absolute value of D^^ \g) is larger than D\^ \g) 
is because, for spherical traps, the dipolar force with neg- 
ative D is more "attractive" than that with positive D. 
In general, the stability boundaries should depend on 
the trap aspect ratio A. However, we find that, for pos- 
itive (negative) D in prolate (oblate) trap, the critical 
dipolar interaction strength is roughly independent of the 
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trap geometry, i.e., D\^^{g) ^ D[ \g) and D\^^{g) ^ 

D[ {g). In fact, under these conditions, the condensate 
shrinks as the dipolar interaction strength grows, such 
that it only occupies a small spatial volume near the sta- 
bility boundary. Therefore, the critical dipolar interac- 
tion strength becomes insensitive to A. For this reason, 
we shall focus on the upper stability boundary {D)^ ^) for 



i(-) 



) for 



oblate traps and the lower stability boundary {D\ 
prolate traps. 

In Fig. [21 we summarize the main results for the stabil- 
ity diagrams in oblate traps. Density oscillations occur 
when the interaction parameters fall into the shaded ar- 
eas, namely, the condensate has a biconcave density pro- 
file in oblate trap. The origin of the density oscillation 
can be intuitively understood as follows: the long-ranged 
dipolar interaction is repulsive in x?/-plane, therefore, 
when the dipolar force exceeds certain threshold value 
the density at the center of the trap is depleted to lower 
the total energy. On the other hand, as one increases 
the dipolar interaction strength D^ the attractive dipo- 
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FIG. 5: (a) Energies as functions of L) for A = 6 and ^f = 25. 
(b) Energies as functions of ^ for A = 6 and D — 55. The solid, 
dashed, and dash-dotted lines correspond to, respectively, the 
kinetic, contact interaction, and dipolar interaction energies. 
The shaded area marks the range of the interaction parame- 
ters where the density oscillation appears. 



lar interaction along 2:-axis also becomes stronger, which 
destablizes the system. In a slightly oblate trap, the sys- 
tem becomes unstable before such threshold dipolar force 
required to from the density oscillation is reached. There- 
fore, for A < 3.2, we do not find any density oscillation 
in condensate wave function for g up to 600. To get a 
better understanding of the density oscillation, we plot, 
in Fig. [3] (a), the effective potential of the system 



t/eff(r) = /7ho(r)+^|^(r)p 



D I dv'Vd{v 
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Apparently, when density oscillation occurs in pancake- 
shaped traps, the effective potential takes a shape of 
Mexican hat. 

The first density oscillations island (denoted as DOI-1) 
appears when A = 3.2. The shape of the DOI clearly indi- 
cates that the stability of the system is slightly increased 
by forming density oscillation in condensate, and away 
from DOI, the D\^ roughly linearly depends on ^. As 
we increase A, DOI-1 moves along the stability boundary 
to negative direction of the ^-axis. In addition, the area 
of DOI-1 shrinks. For even larger A, the second (DOI- 
2) and third (DOI-3) density oscillation regions appear. 
They also move along the stability boundary by following 
the similar pattern as DOI-1 when the trap aspect ration 
is increased. 

A DOI can be roughly characterized using the contact 
interaction strength: the range of ^ parameter covered by 
the DOI, A^DOi, represents its width; while the g value 
at the center of this range, ^doi, is the position of the 
DOI. Figure m shows the trap aspect ratio dependence of 
the positions of the three DOIs; the corresponding width 
is plotted as error bar. One immediately notes that the 
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FIG. 6: The contact interaction strength dependence of the 
upper critical dipolar interaction i^^ for a condensate in 
cigar-shaped traps. The shaded areas mark the DOIs. The 
inset shows the typical result for the isodensity surface of the 
condensate. The broken arrows denote the direction to which 
the DOIs move as A is decreased. 



A dependence of ^doi can be fitted using equation 



^DOi(A) = ae 



-/3A 



7, 



(7) 



where the values of the fitting parameters a, /3, and 7 for 
different DOIs are listed in Tab.|Tl As shown in Fig. 21 the 
fitted curves agree very well with the numerical results for 
all DOIs. Moreover, the values of the fitting parameter 7 
also suggest that increasing A will eventually pushes all 
three DOIs to the region with ^ < 0. For large A, all three 
DOIs are squeezed into a small region of ^ value, but they 
remain disconnected. We point out that even though 
Eq. (jTj) cannot explain why the DOIs appear under those 
contact interaction strengths, it can be used to predict 
the position of DOIs. 

Once the position and width of a DOI is known, one 
may access the DOI by simply increasing the strength 
of dipolar interaction strength to the value close to the 
critical value. In Fig. [5] (a), we present the typical D de- 
pendences of the energies per particle. One immediately 
notices that the energies vary smoothly as the system en- 
ters the DOI, indicating that condensates with and with- 
out density oscillation cannot be distinguished from each 
other by measuring the release energy. As shown in Fig. [5] 
(b), for fixed D^ one may also enters the DOI by lowering 
the contact interaction strength. Again, the energies are 
smooth functions of g across the boundary of a DOI. 

The stability diagram of a dipolar condensate in cigar- 
shaped traps is presented in Fig. [6l where the DOIs ap- 
pear when A < 0.3. In contrast to I) > case, the 
dipolar interaction becomes repulsive (attractive) along 
axial (radial) direction for D < 0. Therefore, along z 
direction, the repulsive dipolar force depletes the den- 
sity at the center of the trap such that the isoden- 
sity surface of the condensate with density oscillation 
takes dumbbell-shaped form. Corresponding to dumb- 
bell structure, the effective potential [Fig. [3] (b)] takes 
the form of a double-well potential. As the trapping po- 
tential becomes more pancake-shaped, the DOIs move 
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FIG. 7: The column density of the expanded cloud for (jj±t = 
(a), 3 (b), 5 (c), and 8 (d). The initial state are obtained 
using (X,g,D) = (0.2,35,-18); during the expansion, g is 
tuned to zero. 



along the stability boundary to the negative direction of 
the ^-axis; in addition, the area of the DOI shrinks. 

We remark that our results do not rule out the possibil- 
ity that there might exist more DOIs for larger g values, 
since our numerical calculations have only covered lim- 
ited ranges of the parameters A and g. 



B. Free expansion 

Now we turn to study the expansion of an initially 
trapped dipolar condensates. Compared the contact in- 
teraction, the anisotropic nature of the dipolar inter- 
action makes the expanded cloud behave quite differ- 
ently [22]. Indeed, it has been used experimentally as 
a diagnostic tool for the detection the dipolar effects in 
Bose-Einstein condensate [4, 8, 9]. 

To study the free expansion dynamics, we consider 
an initially trapped condensate with control parameters 
(X^g^D) falling into a DOI. At time t = 0, we switch 
off the trapping potential, the expansion dynamics of the 
condensate is then described by Eq. ([2]) with Uho = 0. 
If the interactions are ignored completely during expan- 
sion, the condensate will expand ballistically such that 
the time-of-flight image represents the momentum dis- 
tribution of the trapped condensate. From the Fourier 
transform of the initial wave function with density oscil- 
lation, one can easily deduce that two side peaks would 
appear in the time-of-flight image, which can be used 
as the signature of the density oscillation in the initial 
condensate. 

To demonstrate this, we consider the expansion of an 
initially cigar-shaped condensate obtained using the con- 
trol parameter (A, ^, D) = (0.2, 35, —18). In stead of tun- 
ing both g and D to zero, here we only switch off the 
contact interaction for t > 0. Figure [3 shows the column 
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FIG. 8: The column density of the expanded cloud for uj±t = 
(a), 3 (b), and 7 (c). The curves at the upper row show 
the corresponding column density on x-axis, i.e., n(x, 0). The 
initial wave function is obtained using the control parameters 
(A,^,D) = (7,-5,33). 



densities of the expanded cloud, i.e., 



i(x,z,t) = / dyn{x^y^ z^t). 



(8) 



Due to the density oscillation, the initial column den- 
sity has two peaks, in analog to that of the condensate 
trapped in a double-well potential. After the condensate 
expands such that these two peaks overlap, the third peak 
appears as a result of the interference. At time uj±t = 5, 
three peaks have roughly the same height. When the 
system evolves continuously, the heights of the two side 
peaks become lower and lower, but they are still visible 
at uj±t = 8. However, if the contact interaction remains 
unchanged during free expansion, the above scenario will 
be spoiled such that the side peaks vanish very quickly 
for a cigar-shaped condensate. 

Fortunately, for a pancake-shaped condensate, the side 
peaks of the column density remain visible for a long 
time if initially the dipolar interaction is very close to 
the stability boundary. In Fig. [H we demonstrate the 
free expansion of a pancake-shaped condensate with ini- 
tial wave function obtained using the control parameters 
(X^g^D) = (7,-5,33). To make it easier for visualiz- 
ing the structure of the expanded cloud, we also plot the 
column density on x-axis. It can be seen that, even for 
uj±t = 7, the side peaks are still distinguishable from the 
column density, indicating that the time-of-flight image 
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FIG. 9: The upper solid line represents the upper critical 
dipolar interaction strength for a dipolar condensate confined 
in a box potential in the x-y plane and a harmonic poten- 
tial in the z axis. The region between the lower and the 
upper solid lines represents regions where density oscillation 
will occur in the ground state. The inset shows an isodensity 
surface of a dipolar condensate with density oscillation. The 
length L along the x- and ^-axe is taken to be 28a2, where 
ttz = y^K/JrmJz)- The dimensionless interaction strengths are 
defined as ^ = AnNa/az and D = dP / (hw za^z) ■ 



may be used to detect the density oscillation in a dipolar 
condensate. 



V. DIPOLAR CONDENSATE IN A BOX 
POTENTIAL 

In this section we briefly discuss the density oscilla- 
tions of a dipolar condensate trapped in a box potential. 
Specifically, we consider a potential which is harmonic 
along the z-axis with frequency cj^, but is represented by 
infinite square wells with length L along the x- and y- 
axis. Such potential can be realized using, for example, 
tightly focused light sheet |26j . 

Figure [9] displays the stability diagram of such a sys- 
tem. As in the previous case, there exists an critical 
dipolar interaction strength beyond which the system will 
collapse. However, a rather significant difference here is 
the much enlarged region of density oscillation. In fact, 
here there is no disconnected DOIs and density oscilla- 
tion dominates a major portion of the parameter space. 
The inset of Fig. [9] exhibits a typical density plot of a 
condensate trapped in such a potential. One can easily 
see the pronounced density oscillations near the walls. 

That a hard wall induces density oscillations in dipolar 
condensate can be understood from the roton excitation 
of the system [l5[ . Here the effect of the wah is not un- 
like a vortex line [1^ or a localized repulsive potential 
[15|: They all create a local superfluid density depres- 
sion, in the vicinity of which the roton mode becomes 



soft and manifests itself as density oscillations [IS] . The 
same physics also underlies the predicted density ripples 
in superfiuid helium near a boundary [23 or a vortex line 

From this study, one can conclude that density oscil- 
lations in dipolar condensates can be very easily induced 
by a rigid wall. However, in case the condensate is not 
confined in the x-y plane (which can be simulated us- 
ing the periodic boundary condition in the plane), we do 
not observe any density oscillations — as dipolar inter- 
action strength increases, the system remains uniform in 
the transverse plane until it collapses. 



VI. CONCLUSION 

To conclude, we have systematically investigated the 
density oscillations in trapped dipolar condensates by 
mapping out the stability diagram of the system on g- 



D parameter plane. In addition to the biconcave den- 
sity profile in an oblate trap, the dipolar interaction can 
also generate dumbbell-shaped density profiles in pro- 
late traps. On stability diagram, the density oscillation 
occurs only when the interaction strengths fall into the 
DOIs. The relation between the positions of DOIs and 
the trap geometry is explored by fitting the numerical 
results. We also study the free expansion of the conden- 
sate with density oscillation, it is shown that the density 
oscillation can be identified in the time-of-flight image. 
Finally, to improve the numerical precision for the cal- 
culation of dipolar interaction energy, we introduce the 
truncated dipolar interaction potential which is nonzero 
only within a properly chosen cylinder. 
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